#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <time.h>
#include <pthread.h>
#include <stdlib.h>

/* 10.000.000 de dígitos, equivale em bits aproximadamente 33219281. */
#define bits_precisao 33219281

void *calcula_a(void *param);
void *calcula_b(void *param);
void *calcula_p(void *param);
void *calcula_precisao(void *param);

/* Variáveis globais para comparar o erro. */
mpf_t erro_atual, erro_maximo;

/* Variavel para armazenar o pi aproximado com 16 milhoes de casas decimais. */
mpf_t pi_aproximado;

/* Variaveis usadas para o método de Borwein. */
mpf_t a0, a1, a2;
mpf_t b0, b1, b2;
mpf_t p0, p1;

/* Variavel auxiliar para ajudar na troca de contexto entre as iterações e 
	em alguns calculos. */
mpf_t aux_global;

/* Funcao utilizacao para pegar o tempo atual do relogio. */
double tempo()
{
	struct timeval tv;
	gettimeofday(&tv,0);
	return tv.tv_sec + tv.tv_usec/1e6;
}

/* Numero da iteracao atual. */
int iteracao = 0;

/* Variavel que indica de a precisao foi atengida. */
int precisao = 0;

int main(void){

	FILE *output_pi_calculado, *input_pi_aproximado;

	/* Variáveis para registrar o tempo. */
	double tempo_inicio, tempo_fim;

	/* Declaracao das 4 threads. */
	pthread_t thread_a, thread_b, thread_p, thread_precisao;

	/* Seto a precisão padrao para as variaveis a serem inicializadas. */
	mpf_set_default_prec( bits_precisao );


	mpf_init(a0); mpf_init(a1); mpf_init(a2);
	mpf_init(b0); mpf_init(b1); mpf_init(b2);
	mpf_init(p0); mpf_init(p1); mpf_init(aux_global);
	mpf_init(erro_atual); mpf_init(erro_maximo);
	mpf_init(pi_aproximado);

	/* Carrega o valor de pi calculado, do arquivo, na variavel pi_aproximado. */
	input_pi_aproximado = fopen("pi_aproximado.txt", "r");	
	mpf_inp_str(pi_aproximado, input_pi_aproximado, 0);
	fclose(input_pi_aproximado);

	/* Valor do erro máximo, ou seja, número de dígitos corretos. */
	mpf_set_str(erro_maximo, "1e-10000000", 0);

	/* Registra o tempo inicial. */
	tempo_inicio = tempo();

	/* Setando a0. */
	mpf_sqrt_ui(a0, 2.0);

	/* Setando b0. */
	mpf_set_ui(b0, 0.0);

	/* Setando p0. */
	mpf_set_ui(p0, 2.0);
	mpf_add(p0, p0, a0);

	/* Calculando a1. */
	mpf_sqrt(aux_global, a0);
	mpf_ui_div(a1, 1.0, aux_global);
	mpf_add(a1, a1, aux_global);
	mpf_div_ui(a1, a1, 2.0);


	/* Calculando b1. */
	mpf_add_ui(b1, b0, 1.0);
	mpf_mul(b1, b1, aux_global);
	mpf_add(aux_global, a0, b0);
	mpf_div(b1, b1, aux_global);
	
	/* Loop principal das iteracoes. */
	while(1){
	
		/* Cria as 4 thread para executarem em paralelo. */
		pthread_create (&thread_precisao, NULL, calcula_precisao, NULL);
		pthread_create (&thread_a, NULL, calcula_a, NULL);
		pthread_create (&thread_b, NULL, calcula_b, NULL);
		pthread_create (&thread_p, NULL, calcula_p, NULL);
		
		/* Espero a thread que verifica precisao terminal e vejo se ja
			alcançou a precisao desejada. */
		pthread_join(thread_precisao, NULL);
		if(precisao)
			break;
		pthread_join(thread_p, NULL);
		pthread_join(thread_a, NULL);
		pthread_join(thread_b, NULL);


		/* Troco o contexto para preparar para proxima iteracao. */
		*aux_global = *a0;
		*a0 = *a1;
		*a1 = *a2;
		*a2 = *aux_global;
		
		*aux_global = *b0;
		*b0 = *b1;
		*b1 = *b2;
		*b2 = *aux_global;
		
		*aux_global = *p0;
		*p0 = *p1;
		*p1 = *aux_global;

		/*Incrementa o numero da iteracao. */
		iteracao++;
	}

	/* Registra o tempo final. */
	tempo_fim = tempo();

	/* Libero as variaveis alocadas. */
	//mpf_clears(a0, a1, a2, NULL);



	/* Arquivo onde será salvo o resultado. */
	output_pi_calculado = fopen("borwein_paralelo.txt", "a+");

	
	printf("Borwein Paralelo executado com sucesso\n");
	
	/* Salva informacoes sobre a execucao em um arquivo. */
	fprintf(output_pi_calculado, "\n--------------------------------------------------------\n");
	fprintf(output_pi_calculado, " Borwein Paralelo executado com sucesso.\n");
	fprintf(output_pi_calculado, " Os primeiros 30 digitos do resultado sao : ");
	mpf_out_str(output_pi_calculado, 10, 30, p0);
	fprintf(output_pi_calculado, "\n O tempo gasto com a execucao foi de %lf segundos.\n", tempo_fim - tempo_inicio);
	fprintf(output_pi_calculado, " O numero de iteracoes necessarias foi %d.\n", iteracao);

	fclose(output_pi_calculado);


	return 0;
}

/* Funcao da Thread responsavel por calcular An+2. */
void *calcula_a(void *param){

	mpf_t aux;
	mpf_init(aux);

	mpf_sqrt(aux, a1);
	mpf_ui_div(a2, 1.0, aux);
	mpf_add(a2, a2, aux);
	mpf_div_ui(a2, a2, 2.0);

	mpf_clear(aux);

}

/* Funcao da Thread responsavel por calcular Bn+2. */
void *calcula_b(void *param){

	mpf_t aux;
	mpf_init(aux);

	mpf_sqrt(aux, a1);

	mpf_add_ui(b2, b1, 1.0);
	mpf_mul(b2, b2, aux);
	mpf_add(aux, a1, b1);
	mpf_div(b2, b2, aux);

	mpf_clear(aux);

}

/* Funcao da Thread responsavel por calcular Pn+1. */
void *calcula_p(void *param){

	mpf_t aux;
	mpf_init(aux);

	mpf_add_ui(p1, a1, 1.0);
	mpf_mul(p1, p1, p0);
	mpf_mul(p1, p1, b1);
	mpf_add_ui(aux, b1, 1.0);
	mpf_div(p1, p1, aux);

	mpf_clear(aux);
}

/* Funcao da Thread responsavel por verificar a precisao de Pn. */
void *calcula_precisao(void *param){

	mpf_sub(erro_atual, p0 , pi_aproximado);
	mpf_abs(erro_atual, erro_atual);
	if(mpf_cmp(erro_atual, erro_maximo) < 0){		
		precisao = 1;
	}

}
